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Abstract 

Several celestial bodies in co-orbital configurations exist in the solar system. How¬ 
ever, co-orbital exoplanets have not yet been discovered. This lack may result from a 
degeneracy between the signal induced by co-orbital planets and other orbital configu¬ 
rations. Here we determine a criterion for the detectability of quasi-circular co-orbital 
planets and develop a demodulation method to bring out their signature from the ob¬ 
servational data. We show that the precision required to identify a pair of co-orbital 
planets depends only on the libration amplitude and on the planet’s mass ratio. We 
apply our method to synthetic radial velocity data, and show that for tadpole orbits 
we are able to determine the inclination of the system to the line of sight. Our method 
is also valid for planets detected through the transit and astrometry techniques. 


1 Introduction 

Lagrange (1772) found an equilibrium configuration for the three-body problem where 
the bodies are located at the vertices of an equilateral triangle. For relatively small 
eccentricities, the libration around the stable Lagrangian equilibrium points L 4 and 
L 5 is one of the two possible configurations of a stable co-orbital system, called a 
tadpole orbit (by analogy with the restricted three-body problem, we define L 4 as the 
equilibrium point when the less massive planet is 60° ahead of the more massive one 
and L 5 when it is behind). The first object of this kind was observed by Wolf (1906), 
the asteroid Achilles, which shares its orbit with Jupiter around L 4 . At present, more 
than 6000 bodies in tadpole orbits are known in the solar system (MFC 2014). For 
objects in the second configuration, called a horseshoe orbit after the shape the trajec¬ 
tories of the bodies in the corotating frame, the libration encompasses the equilibrium 
points L 4 , L 5 , and L 3 . A single example is known, for a pair of satellites of Saturn 
(see Dermott & Murray 1981b). 

The Lagrangian equilibria points are stable if the masses of the planets are low 
enough. In the quasi-circular case, Gascheau (1843) showed that there is a stability 
condition for the Lagrangian equilibrium 

toqTOi -f TO1TO2 -b mom 2 ^ ^ n nav n'l 

(mo-I-mi-I-TO2)^ 27 

where mo is the mass of the star, and mi and m 2 the mass of the co-orbital planets. 

The mass repartition between the two co-orbitals has a small impact on the stabil¬ 
ity. Within this limit, Gascheau’s criterion guarantees the stability of the linearized 
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equations in the vicinity of L 4 or L 5 . The lower the masses of the co-orbitals with 
respect to the total mass, the larger is the possible libration amplitude. The horseshoe 
domain is stable when the planets have a Saturn-mass or less (Laughlin & Chambers 
2002 ). 

For eccentric orbits, the range of stable mass ratios between the co-orbital and 
the central body decreases as the eccentricity increases (Roberts 2000; Nauenberg 
2002). Moreover, an additional co-orbital configuration exists in the eccentric case, 
called quasi-satellite, as the co-orbitals seem to gravitate around each other in the 
rotating frame. For high eccentricities, co-orbitals have a much larger stable domain 
for quasi-satellites than for tadpole or horseshoe configurations (Giuppone et al. 2010). 

Since the discovery of the first exoplanets (Wolszczan & Frail 1992), a great di¬ 
versity of systems has been found, some of them in mean motion resonances (MMR). 
A few of these resonant systems are highly populated (like the 2/1 MMR), but so 
far no system has been identified in a co-orbital configuration (1/1 MMR). However, 
many theoretical works suggest that co-orbital exoplanets may also exist. Laughlin 
& Chambers (2002) introduced two possible processes that form these systems: (i) 
planet-planet gravitational scattering and (ii) accretion in situ at the L4-L5 points of 
a primary. 

The assumptions made on the gas disc density profile in scenario (i) can either 
lead to systems with a high diversity of mass ratio (Cresswell & Nelson 2008) or to 
equal mass co-orbitals when a density jump is present (Giuppone et al. 2012). In 
their model, Cresswell & Nelson (2008) form co-orbitals in over 30% of the runs. 
Initially in horseshoe configurations, the co-orbital systems are generally damped into 
tadpole configurations. The co-orbitals formed by this process usually have very low 
inclinations and eccentricities (e < 0 . 02 ). 

Lyra et al. (2009) showed that in scenario (ii), up to 5-20 Earth-mass planets 
may form in the tadpole of a Jupiter-mass primary. For existing co-orbitals, the gas 
accretion seems to increase the mass difference between co-orbitals, the more massive 
of the two reaching Jovian mass while the starving one stays below 70Mm (Cresswell 
& Nelson 2009). 

Models that produce co-orbital planets due to dissipation within a disc may also 
experience significant inward migration. Such migration increases the libration ampli¬ 
tude of trojan planets for late migrating stages with low gas friction. This may lead 
to instability, but the damping of the disc usually forces the libration to remain small. 
Equal mass co-orbitals (from super-Earths to Saturns) are heavily disturbed during 
large scale orbital migration (Pierens & Raymond 2014). In some cases, Rodriguez 
et al. (2013) have shown as well that long-lasting tidal evolution may perturb equal 
mass close-in systems. Overall, significantly different mass trojans may thus be more 
common, especially in close-in configurations. 

Co-orbital planets can also be disturbed in the presence of additional planetary 
companions, in particular by a significantly massive planet in another MMR with the 
co-orbitals (Morbidelli et al. 2005; Robutel & Bodossian 2009). Moreover, horseshoe 
orbits are more easily disturbed than tadpole orbits owing to the higher variation of 
frequencies in the system with respect to the libration amplitude (Robutel & Pousse 
2013). 

Detecting co-orbital planets is not an easy task because the signatures of each 
body are usually superimposed. The transit method can solve this degeneracy, either 
by observing both co-orbitals transiting (Janson 2013), or when only one co-orbital is 
transiting by coupling with another detection method (Ford & Gaudi 2006). When 
the libration amplitude is large enough, the transit timing variation (TTV) method 
can also find co-orbital candidates even when only one body is transiting (Janson 
2013; Vokrouhlicky & Nesvorny 2014). 

If the bodies are exactly at the Lagrangian equilibrium point (with no libration), 
the radial velocity (RV) signal is the same as for a single planet, overestimating the 
mass of the system by ~ 13% (Dobrovolskis 2013). Laughlin & Chambers (2002) 
showed that the libration of co-orbitals modulates the RV signal from the star, allowing 
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them to determine a co-orbital system from a simulated RV signal. Cresswell & Nelson 
(2009) estimated that for the co-orbitals obtained in their simulation, this modulation 
is 10% — 20% of the total amplitude of the signal. Giuppone et al. (2012) showed 
that a short-term RV signal (duration inferior to the period of libration) did not allow 
co-orbitals to be distinguished from an eccentric planet or from two planets in 2/1 
MMR (Gozdziewski & Konacki 2006; Anglada-Escude et al. 2010). 

Theoretical and numerical studies seem to agree that tadpole and horseshoe co¬ 
orbitals tend to have low eccentricities and mutual inclinations. In addition, in the 
solar system we observe among the moons of Saturn (Dermott & Murray 1981b, and 
references therein) three co-orbital systems, all with inclinations < 1° and eccentrici¬ 
ties < 0.01. Jupiter’s trojans have a mean inclination of about « 12° and eccentricities 
bellow 0.3. We thus conclude that quasi-circular orbits are a good approximation for 
many co-orbital systems. As the study of quasi-circular and coplanar orbits is easier 
than the full case, in this paper we restrict our analysis to this simpler situation. We 
first give a brief overview of the co-orbital dynamics and an analytical approxima¬ 
tion of the co-orbital motion valid for small eccentricities. In section 3 we derive an 
analytical method that can be used to extract the signature of co-orbitals from the 
motion of the host star. In section 4 we exemplify our method with the radial velocity 
technique using synthetic data, and we conclude in the last section. 

2 Co-orbital dynamics 



Figure 1: Reference angles represented for the circular co-orbital system with respect 
to an inertial frame {x,y). rriQ is the mass of the central star, mi and m 2 the masses of 
the co-orbitals, r* is the distance of the co-orbital i to the central star, and A* its true 
longitude. Following equation (3) one can write A* as a function of Aq, Cj and of the mass 
ratio 5. 


We denote toq the mass of the central star, and mi and m 2 the masses of the 
co-orbital planets. Adapting the theory developed by Erdi (1977) for the restricted 
three-body problem to the case of three massive bodies, and neglecting the quantities 
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of second order and more in 


= 


mi + TO2 
mo + TOi + m2 


( 2 ) 


the mean longitudes Xi and the semi-major axes ai of the co-orbitals can be approxi¬ 
mated by the expressions (see Robutel et al. 2011) 


Ai(t) K. nt Aq , A 2 (t) ~ nt — (1 — d)(^(t) + Ag , 

ai(t) « a (1 - , a2 (t)«a(l-b(l-(5 )|^) , 


(3) 


where n is the averaged mean motion of the barycentre of mi and m 2 , Aq is a constant 
equal to the initial value of the barycentre of the longitudes (Aimi -I-A 2 m 2 )/(mi -|-m 2 ), 
and 


S = 


m2 

mi -I- m2 


(4) 


At first order in /x, we have the constant quantity a = (1 —( 5 )ai-|-i 5 o 2 , which can be seen 
as the mean semi-major axis, and is linked to n by Kepler’s third law = G'(mo -I- 
mi -|- m 2 ), where G is the gravitational constant. Finally, the variable C = ~ -^2 

satisfies the second-order differential equation 
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l -(2 


2cosC)”^/^ 


sinC , 


(5) 


which is one of the most common representations of the co-orbital motion (see Morais 
1999, and references therein). For small eccentricities, it describes the relative motion 
of the two bodies and is valid as long as the co-orbital bodies are not too close to the 
collision (^ = 0 ). 

Equation (5) is invariant under the symmetry () 1 —> — so the study of its phase 
portrait can be reduced to the domain {(,, () S [0, tt] xK (see Fig. 2) (this symmetry will 
be developed later on). The equilibrium point located at {(, () = (tt/S, 0) corresponds 
to one of the two Lagrangian equilateral configurations^, which are linearly stable 
for sufficiently small planetary masses (see below). Another equilibrium point, whose 
coordinates are ( 7 r, 0 ), corresponds to the unstable Eulerian collinear configuration 
of the type L 3 . The separatrices emanating from this last unstable point split the 
phase space in three different regions: two corresponding to the tadpole trajectories 
surrounding one of the two Lagrangian equilibria (in red in Fig. 2), and another one 
corresponding to the horseshoe orbits that surround the three above-mentioned fixed 
points (in blue in Fig. 2). 

As shown in Figure 2, any trajectory given by a solution of equation (5) can be 
entirely determined by the initial conditions {to, Co) such that C(^o) = Co and C(^o) = 0, 
where Co is the minimum value of C along the trajectory, and to the first positive instant 
for which the value Co is reached. 

The possible values of Coj represented by the purple horizontal line in Figure 2, are 
included in the interval ] 0 °, 60°]; Co = 60° corresponds to the equilateral configuration 
where mi is the leading body and m 2 is the trailing one^. The tadpole orbits are 
associated with Co €]Cs,60°[, Cs ~ 23.9° being associated with the separatrix, while 
Co ranges from Cs to 0 for horseshoe orbits. As a result, the shape of the trajectory 
of the relative motion (as the libration amplitude of the resonant angle C) is entirely 
determined by the quantity Co¬ 
in contrast, to and n^/JI are necessary to know the exact position on the trajectory, 
and in particular, the amplitude of the variations of the semi-major axes. We can 

^The coordinates of the other equilateral point are (Stt/S, 0). The permutation of the index 1 and 2 of 
the planets allows the two equilateral conhgurations to be exchanged. 

^The permutation of these two bodies exchanges the two equilateral conhgurations, which are located 
at C = 7r/3 and C = Stt/S, respectively. 
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Figure 2: Phase portrait of equation (5). The separatrix (black curve) splits the phase 
space in two different domains: inside the separatrix the region associated with the tadpole 
orbits (in red) and the horseshoe domain (blue orbits) outside. The phase portrait is 
symmetric with respect to C = 180°. The horizontal purple segment indicates the range 
of variation of Co while the vertical one shows the section used as initial condition to draw 
Fig. 5. See the text for more details. 


rewrite equation (5) by rescaling the time with r = 


cLt^ 


-3 


l-(2-2cosC)”^/^ 


sinC , 


( 6 ) 


where C('r) = C(^)- a consequence, this differential equation does not depend on 
riyfjl. Its solutions are solely determined by the initial conditions tq = y/]Into and 

C('’‘o) = Co = Co¬ 
in a small vicinity of the Lagrangian equilibria, the frequencies of the motion are 
close to 

^0 = 

More generally, excluding the separatrix, the solutions of equation (5) (respectively 
( 6 )) are periodic. The associated frequency, denoted by v (respectively v), depends 
on the considered trajectory. However, the time-normalized frequency associated with 
equation ( 6 ), 

= v/{ny/fl) , (8) 

depends only on Co (11 is plotted versus Co in Fig. 3). In tadpole configurations, this 
dimensionless frequency remains almost constant in the vicinity of the Lagrangian 


5 





















0 10 20 30 40 50 60 


^0 


Figure 3: Variation of the libration frequency u versus Co = Co- The frequency is taken 
over the purple horizontal line in Fig. 2. Inside the tadpole region, the libration frequency 
decreases from y^27/4 at L 4 (Co = 60°) to 0 on the separatrix (Co = Cs ~ 23.9°). In the 
horseshoe domain (Co < Cs) the frequency increases from 0 on the separatrix to inhnity 
when the two planets get closer because the approximations leading to Eq. (5) are not 
valid close to the collision (see Robutel & Pousse 2013). 


equilibrium v Ki vq (Eq. 7) and tends to 0 as the separatrix is reached at Co = Cs- 
In horseshoe configurations, v can take any value. In Figure 3, one can see that far 
from the separatrix, v is always of order unity. This imposes that the variations of 
the difference of the longitudes, Ci are slow with respect to the orbital time scale, i.e. 
V n. It turns out that C,{t)/rL <C 1 and as a consequence, the quantities aj can 
be approximated by a (Eq. 3). Thus, in the circular planar case, at order 0 in ly, the 
position of mi and m 2 in the heliocentric system r = (a; + iy) are given by 

ri = , and rs = ae-hi-<5)Ceh»t+Ao) ^ ( 9 ) 

Within the same approximation, we can also write the derivative of previous equation, 
which gives us the heliocentric velocity of the co-orbitals 

ri = , and h = . ( 10 ) 

While searching for co-orbital bodies, the stability of each configuration also needs 
to be taken into account. In order to determine the influence of the planetary masses on 
the global stability of planar co-orbital systems, we show the results of two numerical 
simulations indicating the width of the stable co-orbital region in different directions. 
In Figure 4, we consider two planets orbiting around a star of mass mo = 1 Mq, with 
fixed initial elements oi = 02 = 1 au, ei = 62 = 0.05, and Ai = zui = 0, and we vary 
the initial element A 2 = ^2 = —Co in [0°,60°] and their masses mi = m 2 = /imo/ 2 , 
with /i/2 G [10“®,10“^]. For each set of initial conditions, the system is integrated 
over 5 Myr using the symplectic integrator SABA4 (Laskar & Robutel 2001) with a 
time-step of 0.0101 year. 

Strongly chaotic systems or systems that quit the co-orbital resonance before the 
integration stops are removed from the computation. In this case, in Figure 4 white 
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Figure 4: Stability of co-orbitals as a function of log]^Q(^) and Co- The initial conditions 
are chosen as to = 0 (Aa/o = 0) and Co £ [0°,60°]: purple horizontal line in Figure 2. 
In black is the separatrix between the tadpole and the horseshoe domain. The stability 
criteria of Gascheau (1843), corresponding to vjn = 11^2, has been indicated. We also 
show the vicinity of two of the main resonances between u and n : the 1/2 resonance (see 
Roberts 2000) and the 1/3. The colour code indicates the value of the libration frequency, 
i.e. logiQ{h'/n). 


dots are assigned to their initial parameters (Co,m)- This strong short-term insta¬ 
bility is mainly due to the overlapping of low-order secondary resonances (Paez & 
Efthymiopoulos 2015). After the elimination of these initial conditions, long-term 
diffusion along secondary resonances may also destabilize the co-orbital systems on a 
much longer time scale. Measuring the temporal variation of the libration frequency 
identihes this diffusion (Laskar 1990, 1999). The black dots indicate a relative varia¬ 
tion of over 10“® between the first and second half of the 5 million years integration 
(to compare with « 10“^° for the long-term stable configurations). They are mainly 
located in the vicinity of the separatrix and near the ejection boundary. In the re¬ 
maining regions, the small variation of the frequency v guarantees, in most cases, the 
stability for a billion years (Robutel & Gabern 2006). For long-term stable systems, a 
colour depending on its libration frequency v is assigned to regular resonant co-orbital 
systems (see the colour code at the bottom of Figure 4). We observe that for large 
planetary masses, slightly lower than the limit value g « 0.037 (Gascheau 1843), the 
stability region is extremely small and strongly perturbed by low-order secondary reso¬ 
nances. The chaos generated by the main secondary resonances, namely the v = n/2, 
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Figure 5: Stability of co-orbitals as a function of logio{fi) and Aa/a. The initial conditions 
are Co = 7i'/3 and Aa/a G [0, 0.06]: vertical pnrple line in Figure 2. The black line indicates 
the separatrix between the tadpole and the horseshoe domains. The dots follow a curve 
Aa oc delimiting the stability region of the horseshoe domain. The colour code 

indicates the value of the libration frequency, i.e. logiQ{v/n) (see Figure 4 for the scale). 


V = n/3, and v = n/4 , shrink the stability region significantly, reducing it to a 
small region near the equilateral configuration (see Roberts 2002; Nauenberg 2002). 
As fi decreases, the width of the stable tadpole region increases, and the destabiliz¬ 
ing influence of the secondary resonances becomes dominant only on the boundary of 
the stability region (see Paez & Efthymiopoulos 2015; Robutel & Gabern 2006; Erdi 
et al. 2007, for the restricted problem). When ^ « 3 x 10“'* « 2Msaturn/MQ, the 
whole tadpole domain becomes stable, except for a small region around the separatrix 
(Co = Cs ~ 23.9°). On the other side of the separatrix, for Co < Cs, stable horseshoe 
orbits start to appear (see Laughlin & Chambers 2002). For lower planetary masses, 
the size of the horseshoe orbital domain increases as g decreases, to reach the outer 
boundary of the Hill sphere at a distance to the separatrix of the order of (see 
Robutel & Pousse 2013). 

In Figure 5 we show another section of the co-orbital region. Instead of varying 
the angle Coi we change the initial value of the difference of the semi-major axes from 
the equilateral equilibrium L 4 towards the outside of the co-orbital region (vertical 
purple line in Fig. 2). More precisely, the initial conditions of the planetary systems 
are ei = 62 = 0.05, Ai = nji = 0, A 2 = W 2 = 7 r/ 3 , and Oj = 1 — (—Ij-^Aa with Aa G 
[0 : 0.06]. As they do in figure 4, the planetary masses vary as mi = m 2 = /.tm.o/2. 
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with ^/2 £ 10“^]. 

The tadpole domaine lies above the solid black line corresponding to the equation 
Aa = 2 1 / 2 / 1 / 3 y//! (Robutel & Pousse 2013). In this case, contrarily to the Co direction 
where the width of the stable tadpole region is a monotonous function of fi, the extent 
of the stability region reaches a maximum for Aa/a « 0.052 at /i = 3.5 x 10“^ and then 
tends to zero with /i as indicated by the above-mentioned curve. For lower values of 
/r, the size of this region decreases until /r reaches the value for which horseshoe orbits 
begins to be stable. After these critical masses, the two domains shrink together but 
at a different rate. The asymptotical estimates of the tadpole’s width in this direction 
is of the order of (black solid line in Fig. 5), while an estimation for the horseshoe 
region is of the order of (black dashed line in Fig. 5, corresponding to the equation 
Aa = 0.47/i^/^) has been fitted to the lower bound of the stable horseshoe region (see 
Robutel & Pousse 2013, for more details). As a consequence, the stability domain 
of the horseshoe configurations becomes larger than the tadpole domain when the 
planetary masses tend to zero (Dermott & Murray 1981a). 


3 One planet or two co-orbitals? 

In some particular situations, co-orbital planets can be identified independently from 
the orbital libration: when both planets are transiting (Janson 2013) or when we 
combine data from transits with radial velocities (Ford & Gaudi 2006). However, in 
general the detection of co-orbitals requires identifying the effect of the libration in 
the data. Vokrouhlicky & Nesvorny (2014) showed that the TTV of only one of the 
co-orbital planets is enough if the libration is large. Laughlin & Chambers (2002) 
showed that the libration induced by co-orbital can have an important effect on the 
radial velocity of a star, while Giuppone et al. (2012) showed that co-orbitals can be 
mistaken for a single planet if the data span is short with respect to the libration 
period. 

In the previous section we saw that co-orbital planets can be stable for large 
libration amplitudes, depending on the parameter /r (see Figures 2 and 4). However, 
the libration period is always longer than the orbital period of the bodies (see the 
colour code in Figure 4). The faster ( librates, the higher the chances of detecting the 
co-orbital bodies, because this reduces the time span needed to detect the libration. 
We write the period associated to the libration frequency iz. The value of P^, 
decreases when /r and n increase (see equation (7) and Figure 3). Therefore, high mass 
ratios and the proximity to the star maximize the detectability of co-orbitals, although 
high mass ratios also lead to the instability of most of the co-orbital configurations 
(Figure 4). Hereafter we consider that the time span of the observations is always 
longer than P^. 

3.1 Signals induced by co-orbital planets 

Most important observational techniques used to detect exoplanets (transits, radial- 
velocity, astrometry) are indirect, i.e. we do not directly observe the planets, but 
rather their effect on the host star. In order to get an idea of the effect of the libration 
of co-orbital planets on the star, we take two examples of co-orbital configurations (see 
Figure 6 ) with the following initial conditions: Ai = 0°, ai = 02 = 1 AU, ei = 62 = 0, 
mi = 0.8 1 O“'^M 0 (red), and m 2 = 1 . 21 O“'*M 0 (blue). In the left graph, ^0 = 25°, 
leading to a large amplitude tadpole orbit, and in the right graph, Co = 23°, leading 
to a horseshoe orbit. The position of the barycentre of the system composed of the 
two planets is represented in purple. With ^ = 2 x 10 “^ and Co near the separatrix, 
these two examples are at the limit of the stability domain, but give a clear idea of 
what we can expect. 

In Figure 7 we show the projection of the stellar orbit on the cc—axis for these two 
configurations. We observe that the signal induced by the Keplerian motion of the 
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Figure 6: Motion of the two co-orbital bodies (red and blue) and their barycentre (purple) 
in a co-rotating frame with frequency n. Tadpole (left) and horseshoe (right). 5 = 0.6. 
Here /r = 2 10“^ and the planets are located at 1 AU from the star. By eliminating the 
influence of n, one can see the long-term motion of the barycentre of the planets. is 
the period of the periodic trajectories represented by the coloured lines. See the text for 
more details. 

co-orbitals is indeed modulated over a period of libration of the resonant angle c). This 
phenomenon was described by Laughlin & Chambers (2002) in the case of a radial 
velocity signal. It is due to the oscillation, with a frequency v, of the distance between 
the barycentre of the two planets and the star, clearly visible in Figure 6. The larger 
the amplitude of variation of the larger the amplitude of modulation. For a given 
Co value, the maximum oscillation amplitude is achieved when 6 = that is, for 
mi = m 2 . In the horseshoe configuration, 5 = 1I2 leads the barycentre of mi and m 2 
to pass by the position of mg, periodically cancelling the signal. 

The bottom panel of Figure 7 shows the spectrum of those signals. The features 
of the spectrum of a modulated signal appear clearly: one peak located at the high 
frequency n and harmonics located on both sides n + pv, where p is an integer and v 
is the frequency of the modulation. In general, the peaks located in n and n ± v are 
the ones with the largest amplitude. However, there is an exception when the signal 
is at the limit of the over-modulation, that is, when the peak located in n disappears. 

This can happen only in the horseshoe configuration when mi « m 2 . In this case, 
the main components of the spectrum would be two peaks separated by 2i^ and the 
system would then be easier to identify. In this paper we focus on the possibility of 
detecting the main three peaks. 

3.2 Motion of a star hosting co-orbital planets 

If the centre of mass of the system is at rest, the position of a star hosting two co-orbital 
planets is given by (in barycentric coordinates) 

ro = ^[(1 - (5)ri-h i5r2] , (11) 

where ri and r 2 are given by equations (9). Since C(t) is a periodic function with 
frequency v, we can expand the terms in Fourier series as 

= E , (12) 

pGZ 
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Figure 7: Motion of the star in the configurations of Fig. 6 in the direction x in the inertial 
frame. In black is the tadpole orbit and in red the horseshoe. The top graph represents 
the evolution of the position of the star over time and the bottom graph its spectrum. In 
these examples, the libration period of the horseshoe orbits is about twice the period of 
the tadpole orbits. See the text for more details. 


where Cp(5, CoTo) is a complex coefficient. Replacing equation (9) and (12) into equa¬ 
tion (11), we get 

ro = , (13) 

pGZ 

with 

Cp = (1 - S) Cp(S, Co, to) + <5 Cp(S - 1, Co, to) , (14) 

and 

ipp = arg{Cp) . (15) 

For the velocity, we thus have (at order 0 in v) 

vq = ifian E \Cp\e^(P''^+^*+^o+Vp) , (-16) 

pGZ 

For instance, if the observational data is acquired through astrometry, we get the 
projection of equation (13) on the plane of the sky, while for radial velocities we use 
the projection of equation (16) in the line of sight. 

The stellar motion can be expressed as the sum of a signal of frequency n, which we 
call the “Keplerian component”, and other signals of frequency n+pv, which we call the 


11 










































































































“modulating components”. For simplicity, we consider only the two main modulation 
components p = ±1, which are the ones with the largest amplitude, hence the ones 
that are easier to detect. We thus introduce the quantity S{t), which represents a 
projection of Tq (equation (11)) or ro (equation (16)) over an observable direction, 
restricted to its main two components, 

S{t)=K{t)+M{t) , (17) 

where 

K{t) = S + So cos(n< + ((>o) , (18) 

and 

M{t) = Si cos{{n + i')t + (j)i) 

+ S'-! cos((n - i/)t + ())_i) , (19) 

where (po, Pi, and p-i depend on and the direction of the projection. Our 

purpose is to check if the Keplerian signal that we have detected is modulated, and if 
our data can be approximated by a signal under the form S{t). 


3.3 Demodulation 

We assume that the Keplerian part of the signal is well determined {Sq and S terms in 
equation (18)), otherwise it would be impossible to look for something else. However, 
the modulating signal {Si and S-i) can be hidden in the noise. In order to isolate the 
effect of the modulation, we suggest using a frequency mixing method similar to the one 
used in the demodulation of radio signals. This method is called “superheterodyne” 
and was introduced by Armstrong (1914). It consists in multiplying the modulated 
signal by a signal that has the same frequency as the carrier. As a result, we obtain a 
peak at the modulating signal’s period in the spectrum. We propose using this method 
on data from co-orbital systems, but it can also be used on any other modulated signal 
produced by a different source (e.g. Morais & Correia 2008). 

We consider a set of N observational data measurements. We denote tk the time 
of each observation and Sk the corresponding observed measurement. First, we fit 
the data with a simple sinusoidal function that contains only the Keplerian part K (t) 
(Eq. 18). This provides us with an initial approximation for S, Sq, n, and pQ. Then, 
we perform a transformation on the raw data Sk to subtract the Keplerian part, 

■Sfe = Sk — K{tk) , (20) 

and then, to isolate the modulation frequency, 

Sk = 4 cos {ntk + (p) , (21) 


where p is an arbitrary phase angle. This modified data set can be fitted with a 
similarly modified function 


where 


S{t) = [S{t) — K{t)]cos{nt + (p) = M{t) cos{nt + (p) 

s 

= cos((2n -I-i/)t-I-(/) 1 -I-(/)) 

s 

+ cos{{2n - I')t + (p-i + (p) 

+ cos(i/t -I- A(p) + AS cos(i/t + (p — (p-i) , 


Si 


= Si cos {(p — (p) 
- (p-i + (pi 




( 22 ) 

(23) 

(24) 
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The libration frequency v is now clearly separated from the Keplerian frequency n. 
As we will see in the following sections, we have A5' <?; The libration contribution 
can therefore be fitted by the term in Si, and the signal is maximized if we are able 
to choose (j) = (j). However, (/> is a priori unknown, so we propose computing the Sk 
for two values of (j) dephased by 7 r/ 2 , for example (p = (po and p = pQ -\- t:/2. By 
proceeding in this way, in the worst case we get p — p = 7 r/ 4 , corresponding to a 
minimum amplitude of Si/\/2. Moreover, by taking the ratio of the fitted amplitudes 
with the two p values, we can additionally estimate p, and thus p±i = pT ^P- 

The initial determination of n using equation (18) always has an error e„, which 
leads to the splitting of the libration term in v into two terms in iz ± e„. Since these 
two frequencies are very close to each other, the Fast Fourier Transform (FFT) usually 
shows a widened peak around preventing an optimal determination of v, Si, and 
Ap. Therefore, once we have some estimations for these parameters, in the last step 
of the demodulation process, we return to the original data set Sk, and directly fit it 
with the full equation S{t) (equation (17)), using the previously determined S, Sq, Si, 
S-i = Si, n, V, pH, pi, and p-i as initial values for the fit. 


4 Detection using the radial-velocity technique 

In this section we apply the general method described previously to the case where 
the data is acquired thorough the radial-velocity technique. In this case, the data 
corresponds to the projection of equation (16) in the line of sight, given by an arbitrary 
direction e*® sin / in the space (Murray & Correia 2011) 

Vr{t) — To ■ e*® sin / = Of ^ \Cp\ cos{pvt + nt-\- pp), (25) 

pGZ 

where 

a = pan sin I , and pp = Pp + 7r/2 — 0 -|- Aq . (26) 

We note that equation (25) could also be the projection of equation (13) over a direc¬ 
tion in the plane of the sky (for example in the case of an astrometric measurement). 
Within our approximation that would only change the value of the parameter a. How¬ 
ever, most of our results on the detectability do not depend on this parameter, thus 
hold true for any measurement technique. For reasons of clarity, we return to the 
example of the radial velocity measurements. 

Considering only the first harmonics of equation (25), one can identify the RV 
signal to the equation (17), which is 

Vr-{t) = S + So cos{nt + po) + S-iCOs{{n — i')t + p-i) 

+SiCOs{{n +iy)t + pi) , (27) 

with Sp = alCpI. We can therefore apply the demodulation process from section 3.3 
to extract the orbital information from the observational data. Our aim now is to 
determine which configurations can be detected for a given precision in the RV ob¬ 
servations, and explain how to retrieve the orbital parameters from the Sp and pp 
parameters. 


4.1 Detectability 

We introduce the following quantity 

, Si + S-I ICil + IC-il 

2^0 2 |Co| 


(28) 


which represents the power of the modulation terms with respect to the Keplerian 
term. When we search for co-orbital planets, the product SoAm must be distinguish¬ 
able from the noise. 
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4.1.1 Detection near the Lagrangian equilibrium 

We consider a system in a tadpole configuration with a low libration amplitude. In 
this case we can use a linear approximation for C near the Lagrangian equilibrium. 
Within this approximation, we can obtain an explicit expression for Vr(t) in terms of 
the orbital parameters. We introduce the small parameter z = — tt/S and write 

C{t) = ^+zcos{iy{t-to)). (29) 

At first order in z and using equation (10), the derivative of equation (11) becomes 


ro = —z/ion 


(1 — d) (l + iSz cos {iz{t — to))) 


+i5(l + z(l — d)zcos(zz(t — to)))e 


„z(nt+Ao+<5f) 


(30) 


Following equation (25), we project equation (30) in the line of sight, and identify the 
terms appearing in equation (27) as 


Sq = aV^ - (5(1 - (5) , 

(31) 

^ ^ ^73(5(1-5) 

ii = i>_i = a----z , 

(32) 

which allow us to compute Am as well: 


. V3 5(1-(5) 

- 2 Vl-5(l-5) ■ 

(33) 


When Too ^ TO 2 > mi the modulation terms can be simplified as 


Si = S-i 


Vs (3 

- zansmi , 

2 mo 


(34) 


where /3 = mim 2 /{mi+m 2 ) is the reduced mass of the planets’ subsystem. We thus see 
that the power in these terms is proportional to P and to the angular separation from 
the Lagrangian equilibrium z. The detection is therefore maximized for large libration 
amplitudes and planets with large similar masses (mi « m 2 ). Nevertheless, we note 
that for planetary systems with mass ratios very different from one, for instance, 
mi/m 2 <C 1, the reduced mass converges to the mass of the smaller planet (/? = mi in 
this case), while for equal mass planets it converges to /3 = mi/2. As a consequence, 
although planets with large equal masses are easier to detect than planets with small 
equal masses, a small mass planet is two times easier to detect if it is accompanied by 
a large mass planet rather than another small similar mass planet. 


4.1.2 Detection in any tadpole or horseshoe configuration 

For large libration amplitudes, we cannot have an explicit expression for Sp with re¬ 
spect to the orbital parameters. Nevertheless, similarly to the linear case, we can prove 
that Am and jCo] depend only on Co and (5. Indeed, since Cp are Fourier coefficients 
of the expression of see equation (12), we can write 

Cp(5, Co, to) = ^ r'"' (35) 

or, in terms of r (see section 2), 

Cp(5,Co,to) = ^ r'’' 
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(36) 








where tq = to/{ny/JI). Since i) depends only on Co, it turns out that 

Cp{S, Co, to) = Cp(5,Co,0)e-*P""« = Cp((5,Co,0)e-*P"*°. (37) 

As a result, the dependence of Cp on tq is explicit. Using the definition of Cp given by 
equation (14), we see that \Cp\, and consequently Am, do not depend on to- 

The dependence of Am and jCol on the parameters (5, Co) is shown in Figures 9 
and 10 for tadpole configurations and in Figure 11 for horseshoe configurations. These 
figures were obtained by integrating the differential equation (5) satisfied by C, with 
initial conditions (C(0),C(0)) = (Co,0). The outputs of these integrations were then 
replaced into the expression of ro for a given set of 5 (Eq. 10). For each simulation, 
the spectrum of a projection of ro has been computed in order to get the value of 
the displayed quantities. These quantities have also been computed from three-body 
direct integrations, which give the same results. 

The RV signal that we obtain for the general cases follows the trends of the linear 
approach. For given values of a, /i, and 6, the detectability of a co-orbital system 
increases as the amplitude of the libration of C increases, i.e. when Co decreases. This 
is still true when Co crosses the separatrix. When S tends to 1 or 0, the modulation peak 
disappears and the signal is similar to the one induced by a single planet. For a given 
Co, Am reaches its maximum when 6 = 1/2. In the horseshoe case, the modulating 
terms have higher amplitudes than the Keplerian term for 0.35 ^ <5 < 0.65, the 
Keplerian term being cancelled when S tends to 1/2. 

We showed at the end of the previous section that a planet of mass rrii (fixed) will 
be easier to detect if its co-orbital companion is significantly more massive (m 2 ^ mi), 
rather than m 2 ~ mi. This holds true in the horseshoe domain, as shown in Appendix 
A.2. 


4.1.3 Detectability for a given data set 

While searching for co-orbital companions of an already detected planet, it is possible 
to put some constraints on what we can expect to observe, based on the observational 
limitations. In addition to the main Keplerian signal, characterized by Ko and P„, we 
also know the time span of the observations, T, and the precision of the instrument, 

e. 

The modulation signal of a co-orbital configuration is detectable if AmKo > e 
(Eq. 28). Thus, the detection of a co-orbital companion can only occur for 


Am ^ 


(38) 


We also know that the libration period is proportional to the orbital period Pn 
(Eq. 8, Fig. 3). One complete libration period can only be contained in the data when 
Pi, > T, therefore 


Pn 


> 


T 


(39) 


The parameter Am depends on S and Co, while the ratio P^/Pn depends on Co and /i. 
The detectability of a co-orbital configuration therefore depends on the mass of both 
planets and on the libration amplitude. 

In Figure 8 we show the ratio Ko/e as a function of the ratio T/P,y, which corre¬ 
spond to the observable quantities. We denote m 2 the most massive of the two planets 
(which is the main contributor to Ko and ^), and mi the mass of the less massive 
planet that we are looking for. We fix m 2 /mo = 10“^ (which is near the maximum 
value allowed for the stability of co-orbital systems) and show the detection limits for 
three different values of mi. Co-orbital companions below each threshold limit can be 
ruled out. 

These detection limits are constrained by the observational limitations (e and P), 
but also by the stability of the co-orbital systems, which is parametrized by the values 
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of Co- A_s Co '^/3 (Lagrange point, with no libration amplitude) or mi/m 2 —>■ 0, 
we have that K^/e > oo. On the other hand, as Co 0, the chances of 

detection increase (the libration amplitude increases), but the system also tends to 
become unstable (Fig. 4). 



Figure 8 : Detectability of a co-orbital companion for 7712/7710 ~ 10 “^. For a given data 
set {Kq! e,TjPn), co-orbital companions with a mass mi can only be detected if they he 
above the respective threshold limit. 


4.2 Characterization of the co-orbital system 

The orbits of the co-orbital planets are fully characterized by the quantities 77, a, 
Co, Ao, to, and sin/. In addition, assuming that the mass of the star is known, we 
can determine the mass of the planets through /r and 5. The frequencies n and v are 
directly obtained when we fit the data with our model (Eq. 27), while a is obtained by 
the third Kepler law from n. Since v depends only on Co (Fig. 3), for each configuration 
there is a bijective map that links ^ and Co given by 




jy 

i'{Co)n ’ 


(40) 


where i> is defined by equation (8). We are thus left with five parameters, Co (or fi), 
S, Ao, to, and sin/, that need to be determined in order to characterize the system. 

We can start looking for the shape of the orbit rather than the exact trajectories 
of the planets as a function of time. Therefore, we ignore by now all quantities that 
depend on Ao and to, i.e. we restrict our analysis to (Eq. 28) and jCol (Eq. 37). 

We define the quantity 47 as 


4' = 2(C) - 4>o) = 4>1+ 4>-i - 2(f)o , 


(41) 
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with (pp = Xq + arg{Cp) + Trj2 — 9 (Eq. 26). Thus 

= arg{Ci{S, Co, to)) + arg{C-i(5Xo,to)) 

-2 arg{Co{S, Co, to)) ■ 

From equation (37), we know that arg{Cp{6Xo,to)) = o,rg(Cp(6, Co, 0)) —pvto- Hence 
'I' depends only on and d. One can show that any quantity defined as a function of 
(pp with p G { — 1, 0,1} and independent of to and Aq is a function of 'k. 

The parameters Am, IOq], and evolve in a different way depending on the orbital 
configuration of the system (tadpole or horseshoe). We thus need to split our analysis 
for these two different configurations. 


4.2.1 Characterization near the Lagrangian equilibrium 

In the linear case, we can entirely determine the trajectories of the co-orbitals analyt¬ 
ically. According to equations (31) and (32), the amplitudes and Si = S-i depend 
on a, Co, and S. By identifying the phases angles appearing in equation (27) to the 
data and then comparing with expression (30), we get three additional equations 


and 


TT \/3S 

())o = Ao-f <5-- arctan(---) , (43) 

6 2 — 0 

(t>±i = K + ^ T t^to . (44) 

3 6 


These three equations, combined with the equations (31) and (32) lead to a system 
of five equations of the form (S'o, 5'i, (/>o, (/'ll <('-i) = Fi<^,S,Co, ^o,to), where F is a 
non-linear function of the five unknown parameters. We can thus get an explicit 
expression for these parameters from the observational data. Then, the expression of 
v near the Lagrangian equilibrium (Eq. 7) can be used to get the value of p. Finally, 
the inclination I can be deduced from the definition of a (Eq. 26) : 


sin / = —^ . (45) 

pan 

We can thus remove the classic p sin / degeneracy in this case and fully determine the 
exact masses of the planets and their trajectories in space. 

Replacing expressions (43) and (44) for pp in the expression of 'k (Eq. 41) gives 


/ . TT 

4- = 2arctan(^^) - - , 


(46) 


i.e. near the Lagrangian equilibrium 4' only depends on 6. Since 0 < d < 1, we have 
41 G [—7r/3,7r/3], and for S = 1/2 we get 4' = 0, which corresponds to equal mass 
planets. 


4.2.2 Large amplitude tadpole orbits 

As discussed in section 4.1.2, for large libration amplitudes it is not possible to obtain 
an explicit expression for the orbital parameters from the Sp terms. The same applies 
to the phase angles pp. However, for tadpole configurations it is still possible to 
inverse the problem using implicit functions and to fully characterize the orbits from 
the modulation terms in equation (27). 

In Figure 9, we show iso-values of the parameters Co and <5 with respect to the 
quantities Am and 4/ (see section 4.1.2 for more details). For tadpole orbits, we see 
that each couple {Am, 4^) corresponds to a unique couple ((o, <5). One can thus 
determine the values of Co and 6 directly from Am and 4^. 
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^y[rad] 



Figure 9: Level curves of 6 (black) and (q (red) for the tadpole configuration, with respect 
to Am and 'k. See the text for more details. 


We also know that |Co(^, Co)l depends only on (q and <5 (see section 4.1.2). In 
Figure 10 we show iso-values of ICo]. Since So = a|Co|, we can directly obtain the 
value of a from (^O) 5), and hence from {A^m 4'). We can thus determine sin J (Eq. 45), 
since /i is linked to (g through expression (40). The parameters S, (o, and sin J are 
then fully determined for the tadpole configuration. 

Finally, similarly to the linear case (section 4.1.1), for a given 6 one can show that 
(pQ is a bijective map for Ag G [0, 27r/n[, and 4>i — Ag is a bijective map for tg S [0, 27: /v[ 
(see equations (43) and (44)). The values of Ag and tg are therefore determined by 
the values of (pg and (pi. Then, one can use equations (3) and (5) to obtain the orbital 
parameters of the co-orbitals. 

4.2.3 Horseshoe orbits 

For the horseshoe configuration, it is also not possible to obtain explicit expressions 
for the orbital parameters from the Sp terms. However, a symmetry in (p allows us to 
compute this (see Appendix A.l): 

T = arg{Ci) + arg{C-i) - 2arg{Cg) = tt . (47) 

Since T is constant in horseshoe configurations, we cannot use it to get an additional 
constraint on the couple (<5, Co)- 

In Figure 11 we plot [Cg] and A^ICgl = (ICil -|- |C'_i|)/2 versus 5 (see section 4.1.2 
for more details). The graphs are symmetric in <5 = 1/2. One can see that these quan¬ 
tities vary significantly with 5, but are are almost constant in regard to Co (different 
colour curves in Fig. 11), except near 5=1/2 for AmlOgj. Thus, we can assume an 
average value for Co in the horseshoe domain. From this average value, we get approx¬ 
imated values of 5 and a by knowing Am and ICg/ Then, we can obtain approximated 
values for the parameters tg and Ag from (pg and (pi^ as explained in the tadpole case. 
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Co 


Figure 10: Level curves of ICqI for the tadpole configuration with respect to Co and S. See 
the text for more details. 




6 6 

Figure 11; Left: ICqI with respect to 6 in the horseshoe configuration. As Am,(S = 1/2) = 
+ 00 , we plot the quantity Am\Co\. Right: Am\Co\ versus 6 in the horseshoe configuration. 
These quantities are symmetric with respect to (5 = 0.5. red: Co = 23°, purple: Co = 19°, 
blue: Co = 15°. See the text for more details. 
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However, the degeneracy in fi sin I remains, because of the strong dependence of v on 
Co in the horseshoe domain (see Fig. 3). One of the ways to get this information is to 
consider higher order harmonics in the expansion of the radial velocity, equation (25). 
However, as these harmonics are about 10 times smaller than Si, much more accurate 
data is required. 

4.2.4 Tadpole or horseshoe? 

Since the method that we use to determine the orbital parameters of a co-orbital sys¬ 
tem depends on its conhguration (tadpole or horseshoe), it is legitimate to ask whether 
it is possible to know the configuration type before we choose one method or another 
for reducing the observational data. 

Once the signature of a co-orbital system is detected (by the observation of a mod¬ 
ulation in the radial-velocity data) we can compute 4* from equation (42). One can 
see from equation (58) that 4' = tt in the horseshoe configuration, while 4^ G [~2j 2] 
in the tadpole configuration (Fig. 9). Since the domains for 4t are exclusive in the 
different configurations, by computing 4' we can immediately distinguish between a 
horseshoe and a tadpole configuration. 

When the detected signal is at the limit of the instrumental precision, the phases 
can be improperly determined. In this case, one can always compute using ex¬ 
pression (28). As shown in Appendix A.2, ranges within [0, -|-oo[ in the horseshoe 
configuration. In the tadpole configuration. Am reaches its maximum value for (5 = 1/2 
and Co near the separatrix. We can see in Fig. 9 that this quantity remains below 1/3. 
Therefore, Am > 1/3 is also a sufficient condition to know that a co-orbital system is 
in a horseshoe configuration. 

4.3 Application to synthetic data 

We now apply the methods developed in the previous sections to two concrete situa¬ 
tions of stars hosting a pair of coobital planets in quasi-circular orbits, one for tadpole 
and another for horseshoe orbits. In Table 1 we list the initial osculating orbital ele¬ 
ments for these two hypothetical systems orbiting a solar-mass star. We then generate 
synthetic radial-velocity data for these systems by numerically integrating the equa¬ 
tions of motion using an n-body model. In order to create a realistic data set, we 
use the same observational dates taken for the HD 10180 system (Lovis et al. 2011) 
to simulate the acquisition days, and associate with each measurement a Gaussian 
error with cr = 1 m/s. These synthetic data sets contain 160 measurements spanning 
4600 days and correspond to an instrumental precision of ~ 1 m/s. The orbital peri¬ 
ods of the planets are around 11.5 days in both examples, such that we can observe 
at least three complete libration cycles over the length of the observations. 

4.3.1 Tadpole orbits 

Our tadpole system is composed of two Saturn-like planets at 0.1 au (comparable 
masses and eccentricities). The individual RV amplitudes of both planets are K ~ 
10 m/s, well above the instrument precision. Therefore, the signatures of the planets 
can be easily identified in the data, and we use this example to illustrate how to 
retrieve the complete set of orbital parameters listed in Table 1 with our method. 

In Figure 12(a), we show a generalized Lomb-Scargle periodogram of the radial 
velocity data (Zechmeister & Kiirster 2009). The Keplerian component of the signal 
with a period « 11 days can clearly be identified. We fit the raw data with a 
Keplerian function (Eq. 18) and obtain an initial estimation for P„ « 11.46 days, 
S « 6.5 m/s. So ~ 61.1 m/s, and ~ 341.6°. We then subtract the Keplerian 
contribution to the data and obtain a modified data set s/ (Eq. 20). In Eigure 12 (b), 
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Table 1: Osculating orbital elements for a given date of two hypothetical co-orbital systems 
orbiting a solar-mass star. 


par am. 

tadpole 

planet 1 planet 2 

horseshoe 

planet 1 planet 2 

m [M®]... 

a [au]. 

A [deg].... 

e. 

w[deg] ... 

7 [deg].... 

S [km.s~^] 

200 100 

0.0987 0.1013 

0 300 

0.05 0.05 

0 300 

60 60 

6.500 

17.15 3.00 

0.1 0.1 

0 339 

0 0 

0 0 

90 90 

6.500 

a [au]. 

Co [deg] ... 

5. 

0.09955 

37.00 

0.3333 

0.10000 

21.00 

0.1488 


we show a periodogram of this modified data. We observe that the main peak with 
a period of approximately 11 days is replaced by two nearby smaller peaks. This is a 
clear indication of the presence of a modulation, each peak corresponding to the n±v 
terms (Eq. 27). 

In order to better determine the libration frequency, we modify the data again 
following expression (21) adopting (j) = (j)Q = 341.6° and (j) = (j)o + 7r/2 = 71.6°. In 
Figures 12 (c) and 12 (d) we show the periodograms of Sk for these two transformations, 
respectively. In both transformations we observe that the peak around 11 days is 
replaced by some power at the periods near 5 and 150 days, corresponding to the 
frequencies 2n and v, respectively (Eq. 22). However, while for (p = (j)Q the maximum 
power is observed for i/ (Fig. 12c), for p = pQ + -kI2 it is observed for 2n (Fig. 12 d). 
From expression (23), we see that the amplitude Si associated with the term with 
frequency v is reduced by 

Si{(j)) = Sicos{p-(p) = Sicos(^ + (pn-. (48) 

For tadpole orbits we have ~ 0 (Fig. 9), which means that p ^ po (Eq. 41). There¬ 
fore, Si is maximized for p ^ pQ and minimized for p ^ pQ+'iTj2 (Eq. 48). Performing a 
FFT to Sk allows us to estimate Pi, « 154.66 days, Si « 4.23 m/s, and Ap « —116.5°. 
We can also estimate p (and hence pi and p-i) using the ratio between the two am¬ 
plitudes 

P = Po+ arctan ^ 

^ SiiPo) ^ 

Finally, adopting these values as initial parameters, we refit the raw data Sk by per¬ 
forming a minimization of expression (27) using the Levenberg-Marquardt method 
(e.g. Press 1992). The results corresponding to the minimum of are shown in 
Table 2. 

From the observational parameters listed in Table 2, we can obtain the corre¬ 
sponding orbital parameters using the inversion method explained in section 4.2.2. 
The osculating orbital elements are then obtained through the equations (3) and (5). 
The results are given Table 3. Except for the eccentricities and the longitudes of the 
pericentre, which cannot be determined with a Keplerian circular orbit approximation 
(Eq. 18), we obtain a very good agreement for the remaining parameters (cf. Table 1). 

We can still improve the quality of the fit in a last step, by performing an ad¬ 
justment to the data using the direct n-body equations of motion (e.g. Correia et al. 
2010). By adopting the orbital parameters listed in Table 3 as the starting point, the 
algorithm converges rapidly to the best ht. The results are given in Table 4. This 
last step slightly improves the orbital parameters obtained previously (lower and 
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Figure 12; Periodograms of the synthetic radial velocity of the tadpole configuration 
presented in Table 1. (a) raw data Sk] (h) modified data after the subtraction of the 
Keplerian signal (Eq. 20); (c) modified data Sk with (p = cpo (Eq. 21); and (d) modified 
data Sk with (p = cpQ + 7r/2 (Eq. 21). 


rms), because it is able to additionally fit the eccentricities and the longitudes of the 
pericentre. We note, however, that the n-body algorithm is only able to converge to 
the correct orbital solution because it used the parameters from Table 3 as starting 
point. Indeed, the phase space of co-orbital planets has many other local minima that 
provide alternative solutions that are not real. 

4.3.2 Horseshoe orbits 

Our horseshoe system is composed of a Neptune-mass and a 3 Earth-mass planet at 
0.1 au. It is at the limit of detection, since the individual RV amplitudes of each planet 
are K = 4.85 m/s and K — 0.85 m/s, respectively. With this example we intend to 
show the limitations of our method. 

In Figure 13 (a), we show a generalized Lomb-Scargle periodogram of the RV data. 
As for the tadpole example in the previous section (Fig. 12), the Keplerian component 
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Table 2: Fitted parameters using expression (27). 


param. 

tadpole 

horseshoe 

Pn [day]... 

11.4599 ±10"“ 

11.5492 ±5 10"'‘ 

[day]... 

154.66 ±0.06 

1340 ±19 

S ]km/s] .. 

6.5001 ±10~'‘ 

6.5001 ±10"'^ 

5o [m/s]... 

61.1±0.1 

4.9±0.1 

5i [m/s]... 

4.23 ±0.09 

1.2±0.1 

5-1 [m/s] . 

4.23 ±0.09 

1.2±0.1 

Po [deg] ... 

341.6 ±0.1 

22.93 ±1.69 

pi Jdeg] ... 

266.4 ±1.8 

309.3 ±6.4 

P-i [deg].. 

33.3±2.1 

280.2±7.9 

yp . 

2.570 

1.613 

rms[m.s“^] 

2.8217 

1.8489 

Am . 

0.069 

0.247 

Vk [deg].... 

-23.5 

183.64 


Table 3: Osculating orbital elements obtained through the inversion of the harmonic 
terms fitted to the observational data (Table 2). The elements marked with * cannot be 
determined with the Keplerian circular orbit approximation (Eq. 18), so they have been 
fixed at constant values. ** indicates that the displayed mass is the lowest possible value 
(msinl). 


param. 

tadpole 

planet 1 planet 2 

horseshoe 

planet 1 planet 2 

m[M(s] 

226.4 

101.6 

18.19** 

2.74** 

a [au].. 

0.099 

0.101 

0.1002 

0.0987 

A [deg]. 

1.380 

303.5 

5.10 

320.78 

e. 

0* 

0* 

0* 

0* 

ro [deg] 

0* 

0* 

0* 

0* 

I [deg]. 

59.85 

59.85 

* 

o 

90* 

a [au].. 

0.09948 

0.09999 

Co [deg] 

38.01 

18.5* 

5. 

0.3440 

0.1309 


of the signal can clearly be identified for a period « 11 days. We thus fit the 
raw data with a Keplerian function (Eq. 18) obtaining an initial estimation for P„ « 
11.55 days, S « 6.5 km/s, Sq « 4.9 m/s, and (/>o ~ 23°, subtract its contribution to 
the data, and obtain a modified data set s/ (Eq. 20). However, unlike the tadpole case, 
in the new periodogram of the residual data, there is no clear peak above the noise 
(Fig. 13 b). Therefore, such a system can easily be mistaken with a system hosting a 
single planet at 11 days. 

We can nevertheless apply our method to search for the traces of a co-orbital 
companion. We thus modify the data Sk according to expression (21) adopting (p = 
po = 23° and p = po + 7r/2 = 113°. In Figures 13(c) and 13(d) we show the 
periodograms corresponding to these transformations, respectively. For p = pQ the 
periodogram is very similar to the one with the residual data (Fig. 13 b), so we conclude 
there is nothing else above the noise in the data. However, for p = (/)o+7r/2 the scenario 
is completely different as a significant peak appears around 1500 days, corresponding 
to the libration frequency (Fig. 13d). Indeed, for horseshoe orbits we have 'k = tt 
(E q. 58), which means that p = po + 7r/2 (Eq. 41). Therefore, Si is null for p = pQ 
and maximized ior p = pQ + ir/2 (Eq. 48). 

Performing a FFT to Sk allow us to estimate P,y « 1340 days, « 1.2 m/s, and 
Ap « 295°. Adopting these values as initial parameters, we refit the raw data Sk 
with expression (27). The results corresponding to the minimum of 3,re shown in 
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Table 4: Best fitted orbital solution using the direct n-body equations of motion, and 
adopting the orbital parameters listed in Table 3 as the starting point. 



tadpole 


horseshoe 

param. 

planet 1 

planet 2 

param. 

planet 1 

planet 2 

m [M®] ... 

195.68 ±0.31 

100.40 ±0.35 

m [M®] ... 

18.79 ±0.008 

2.99 ±0.005 

a [au]. 

0.099 ±610"® 

0.101 ±110“'‘ 

a [au]. 

0.100 ±210“® 

0.099 ±710“® 

A [deg] .... 

2.3T1.7 

306±5 

A [deg] .... 

4.69±2.83 

318.42±8.56 

e . 

0.056T0.003 

0.049±0.003 

e . 

0.000±10“® 

0.000±10“® 

ro [deg].... 

0.01T0.04 

304.2±0.4 

w [deg].... 

0.000±10“® 

0.000±10“® 

7 [deg] .... 

65±2 

57±3 

I [deg] .... 

90* 

90* 

a [au]. 

0.09953 

a [au]. 

0.1000 

Co [deg].... 

37.40 

Co [deg].... 

21.57 

5 . 

0.3391 

5 . 

0.1371 

. 

1.652 

yp . 

1.595 

rms[m.s ^] 

1.8770 

rms[m.s“^] 

1.8513 


Table 2. Comparing these results to the tadpole case, we observe that the uncertainty 
associated with the Sp and (j)p terms is larger, but still near 1 m/s, which corresponds to 
the considered precision of the instrument. Our method is therefore able to extract any 
information on the existence of a co-orbital companion, provided that the information 
on the libration terms is accessible in the data. 

Once the existence of a co-orbital companion is confirmed, we can determine the 
corresponding orbital parameters. The parameter i/o (which gives the departure of the 
semi-major axis and the mean longitudes from their mean value) has a low impact on 
the orbital parameters and cannot be easily determined in horseshoe configuration (see 
section 4.2.3). However, its value is constrained by the stability of the system: in the 
horseshoe configuration, it ranges between its lowest stable value for M sin /, in 

our case « 6 x 10“^ (see Fig. 4), and the separatrix. We therefore have G [13°, 24°]. 
We take Co = 18.5° (average value on this interval) and compute the corresponding 
orbital parameters. We obtain a system close to the original one (Tab. 3). 

In the horseshoe case, we cannot determine either the eccentricities and the lon¬ 
gitudes of the pericentre, because we used a Keplerian circular orbit approximation 
(Eq. 18), or the inclination to the line of sight, because we only fit the first three 
harmonics (Eq. 27). In a final step, we perform an adjustment to the data using the 
direct n-body equations of motion, and we obtain a similar adjustment (Tab. 4). 


5 Discussion and conclusion 

In this paper we have revisited the dynamics of quasi-circular co-orbital planets. By 
computing their gravitational effect on the parent star, we have found a simple method 
for detecting this kind of planets, provided that the orbital libration can be seen in 
the observational data. Indeed, when the star is accompanied by co-orbital planets, 
in addition to the Keplerian orbital motion, there is a modulation at a longer period, 
corresponding to the libration frequency. Therefore, commonly used methods for 
signal demodulation (see section 3.3) can also be applied to co-orbital systems, allowing 
the amplitude and the frequency of the modulation to be identified more accurately. 

Every time a modulation is observed in the motion of a single planet, an inquiry 
should be made to check if it can correspond to the libration induced by another 
co-orbital planet. In this paper, we explain a way to quantify which co-orbital con¬ 
figurations can be expected: for stability reasons, we can put boundaries for pairs of 
the parameters (/i, Co); for data span duration reasons, we can estimate the frequency 
of libration v, depending on (n, /i, Co); for measurement precision reasons, we can 
estimate the amplitude of the modulating peaks, which depends on the parameters 
(/r, n, (5, Co)- 
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Figure 13: Periodograms of the synthetic radial velocity of the horseshoe configuration 
presented in Table 1. (a) raw data s^; (b) modified data after the subtraction the 
Keplerian signal (Eq. 20); (c) modified data Sk with (p = cpo (Eq. 21); and (d) modified 
data Sk with (p = cpQ + 7r/2 (Eq. 21). 


For reasons of clarity, we exemplify our method in the case of a radial velocity 
signal. However, our results are valid for any other method that measures a projection 
of the stellar motion. We have shown that the relative amplitude of the modulation 
signal depends only on the distance to the Lagrangian equilibrium, ^Oj and mass 
ratio, 5. Therefore, the detection of co-orbital planets is enhanced for large libration 
amplitudes around the Lagrangian equilibrium (i.e. small Co values), and for planetary 
masses equally distributed between the two co-orbitals (d « 1/2). 

In order to reduce the data, we proposed a direct inversion from the periodograms 
of the signal to the osculating elements of the system. For systems in the tadpole 
configuration we are able to determine the inclination of the orbital plane with re¬ 
spect to the plane of the sky and hence the true masses of the planets (and not only 
the minimum masses). In the horseshoe case this is not possible without considering 
higher harmonics for the modulation. 
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A Symmetries 

Equation (5), respectively (6), possesses several symmetries. We use two of them to 
study analytically some features of the horseshoe configuration. On the one hand, we 
have the symmetry with respect to (^ = 0 (Aa/a = 0 in Fig. 2): 


C{-{t - to)) = ({t - to) . 


(50) 


On the other hand, we have the central symmetry of the phase space (^, Aa/a) in 
C = TT and Aa/a = 0: 



Similar expressions can be obtained for (/. In the tadpole configuration, these sym¬ 
metries exist as well, but the symmetry (51) maps a vicinity of L4 to a vicinity of 
Lb- 

We can use these symmetries to simplify the expression of the coefficients Cp given 
by equation (14). Our purpose is to study the values of Am{S, Co) and 4'((5, Co)- Since 
none of them depends on the value of tq, we take Tq = 0 from now on. The coefficients 
Cp (Eq. 36) become 



(52) 


Since we took tq = 0, is an even function in the case of a horseshoe orbit. Hence 
g-ipjyr becomes cos{pvT) in the expressions of the Cp. By splitting this expression into 
two integrals and changing t io t + t:/ u m. the first one, we get 



(53) 


Then, using the symmetry given by expression (51), the previous integral simplifies as 



(54) 


hence 


As a consequence, using equation (14), we get for p = 0 



(55) 



(56) 


5cos((5 — l)(7r — C(t))) dr. 


and for q = ±1 



( 57 ) 


5 sin(((5 — l)(7r — C,{t))) cos{qi>T)dT . 


We obtain Ci = C-i. 
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A.l Computation of ^ 

From equation (56), we have arg{Co{S)) = i57r if 5 € [0, l/2[ and + tt if <5 €]l/2,1]. 
Since arg{Ci) = arg{C-i) = {S — l/2)7r (Eq. 57), we conclude that for any horseshoe 
configuration 

= arg{Ci) + arg{C-i) - 2arg{Co) = tt , (58) 

i.e. is constant and equal to tt. 


A.2 Computation of 

Generally, the \Cq\ (Eq. 57) does not have an explicit expression. However, Am can 
be computed for some specihc values of 5. We denote = Cq{6XoA)- For 6 = 1/2, 
we have 


Ci/2 



sin((7r - C(t))/2) 


cos{qvT)dT . 


(59) 


The amplitude of the first harmonics {g = ±1) of the Fourier series of an odd function 
is not null. Thus, since from expression (56) we have that \Cq \ = 0, we can conclude 
that in the horseshoe configuration Am{\Ao) = oo (Eq. 28). Similarly, one can also 
see from equations (56) and (57) that Ar„(0, Co) = ^m(l, Co) = 0. 


B Mass ratios 


In section 4.1.1, we have shown that in the vicinity of the Lagrangian equilibrium, 
a planet with mass mi is easier to identify when its co-orbital companion is much 
more massive (toi m 2 ) rather than when mi « m 2 . We show here that this result 
holds true in the horseshoe configuration. Using the symmetries (50) and (51), one 
can rewrite Ci (Eq. (57)) as 


Cl = — 
TT 


(l-(5) sin(<5(7r-C(— - T-))) 


(60) 


—(5sin(((5 — l)(7r — C(- t))) smlvAdr. 

2v J 


Eor a mass mi, we want to compare the quantity Si = AmSo = a\Ci\ in the case 
of mi = m 2 {6 = 1/2) against the case when mi <§; m 2 (5 « 1 — mi/m 2 = 1 — e). 
Writing X{t) = tt — — r), from equation (57) we have 


\cC = 


2v 

TT 


[sin(X(T)/2)] sin(i>r)(i7 


and at first order in e, equation (57) yields 


|C/-| = 


97 / 

e — / [X{t) sm{X{r))]sm{i>r)dT 


(61) 


(62) 


We have /f(0) =0 and /f(^) = tt — Co- One can show that X is a monotonous 
function in the interval r € [0,77/(2!/)], hence sin(/f/2) and X -|-sin(X) are a positive 
function in this interval. Moreover, for X € [0,7r], we have the following inequality: 


7 rsin(/f/2) < (X -|-sin(X)) < 4sin(M/2) . 


(63) 


Since sin(z/r) is also a positive function on the considered interval, the inequality in 
equation (63) holds true when we multiply each term by sin(i/T) and integrate over 
T e [0,7r/(2i/)]. Einally, we get 

ttC^^ < C\-^/e < 4Cy^ (64) 
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When 6 = 1/2, we have fj, « 2milmQ, while when <5 = 1 — e, we get fi « mi/(emo). 
Multiplying equation (64) by a, we obtain 

< Sl-"^ < 2Sy^. (65) 

We finally conclude that in the horseshoe case, for a given mass mi, the co-orbital 
couple (toi,TO 2 ) is up to two times easier to identify when mi ^ m 2 rather than when 
mi « m 2 . 
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